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We develop a continuum elastic approach to examining the bending mechanics of semiflexible filaments with 
a local internal degree of freedom that couples to the bending modulus. We apply this model to study the 
nonlinear mechanics of a double stranded DNA oligomer (shorter than its thermal persistence length) whose 
free ends are linked by a single standed DNA chain. This construct, studied by Qu et al. [Europhys. Lett., 94, 
18003, 201 1], displays nonlinear strain softening associated with the local melting of the double stranded DNA 
under applied torque and serves as a model system with which to study the nonlinear elasticity of DNA under 
large energy deformations. We show that one can account quantitatively for the observed bending mechanics 
using an augmented worm-like chain model, the helix coil worm-like chain. We also predict that the highly 
bent and partially molten dsDNA should exhibit particularly large end-to-end fluctuations associated with the 
fluctuation of the length of the molten region, and propose appropriate experimental tests. Wc suggest that the 
augmented worm-like chain model discussed here is a useful analytic approach to the nonlinear mechanics of 
DNA or other biopolymer systems. 

PACS numbers: 87.14.gk, 87.10.Pq, 87.15.hp, 87.15.La 



I. INTRODUCTION 

The worm-like chain (WLC) [1] has become the standard 
model of the mechanical behavior of a variety of nanoscale 
filamentous structures in thermal equiUbrium. This approach 
relies on a single elastic parameter, the bending modulus of 
the filament, and posits that the bending response to an ap- 
plied torque is simply linear. The well-known and dramatic 
force vs. extension nonUnearities [2] of the WLC are funda- 
mentally geometric in origin, originating from the inextensi- 
bility of the filament. The bending mechanics of carbon nan- 
otubes, microtubules, DNA, and other filamentous materials 
can all be modeled as simple elastic rods as long as the lo- 
cal bending deformations are sufficiently small and thus their 
behavior in equilibrium can be well described by the WLC. 

More recently, it has become clear that, for sufficiently 
high energy deformations, there are significant discrepancies 
between the WLC predictions and the observed mechanics 
of such filaments. For example, carbon nanotubes (CNTs) 
do not behave as linearly elastic rods at large stresses and 
strains. When bent sufficiently sharply, the competition be- 
tween bending and stretching in multi-walled CNTs induce 
rippling modes of deformation [3-9] as shown in Fig. lb. The 
resulting kinks that form in sharply bent CNTs are reminis- 
cent of the buckled shape that occurs in highly bent drink- 
ing straws or tape-measures, a phenomenon known as Brazier 
buckling flO, 11]. 

Similarly, there is the now significant evidence that 
biopolymer filaments display more complex elastic behavior 
that cannot be understood in terms of a simple WLC model. 
This includes the work on DNA structural transitions induced 
by applied tension [12-19], which have been modeled as a 
type of two-state, or helix coil transition [20]. More germane 
to the current issue are the studies of bending DNA that can 
be explained only by assuming that the bending modulus of 
DNA is anomalously low in higher energy bends, i.e., ones 



having a sufficiently small radius of curvature [21-28]. One 
reasonable interpretation of these data that is consistent with 
the mechanical studies of bent carbon nanotubes and macro- 
scopic filamentous structures is that at high applied stresses 
the system accesses other internal modes of deformation (e.g. 
rippling of CNT), which, by increasing the number of accessi- 
ble degrees of freedom, necessarily increases the compliance 
of the filament towards bending. This leads to a nonlinear, 
stress-induced softening of the material. 

For dsDNA in particular, we recognize there are many such 
internal degrees of freedom associated with twist, roll, and 
slide motions of base pairs [29], but we expect that the princi- 
pal effect leading to a dramatic reduction in the local bending 
modulus of the double stranded helical structure is the break- 
ing of hydrogen bonds associated with the base pairing inter- 
action and the consequent "flip-out" of these individual bases. 
In essence, we expect high curvature regions of the chain to lo- 
cally melt into effectively two weakly coupled single stranded 
chains, which are significantly more compliant to bending. 

It may be possible to make progress on developing more de- 
tailed theories of DNA mechanics at short length scales and in 
high elastic energy configurations by including atomistic de- 
tail in order to capture these effects. This approach, however, 
sacrifices the elegance and great utility of the WLC model, 
which accounts for the statistical mechanics of a variety of 
semiflexible filaments using a minimal set of material param- 
eters - a single bending modulus. We propose that there is 
a useful intermediate level of description, which serves as a 
minimal extension of the WLC that allows for the observed 
elastic nonlinearities. 

Our approach is to treat the filament as a two-state elas- 
tica [22, 23, 30-34] in which the filament has a higher bend- 
ing modulus in the low energy state as compared to the higher 
energy one. For example the melted, more compliant regions 
of double-stranded DNA (dsDNA) have a higher free energy 
per unit length than the ordered dsDNA, but have a signifi- 
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cantly lower bending modulus. This approach to the elasticity 
of DNA, featuring two local internal states coupled to the con- 
formation of the chain, has been explored previously in studies 
of the B-DNA to S-DNA transition under tension [12, 20] and 
in the study of DNA cyclization, where the same local DNA 
melting as we propose here was implicated in the enhance- 
ment of cycUzation probabihties [22]. A similar description 
of the kinked region of a CNT applies. This assumption leads 
generically to the nonlinear softening of the filament at suffi- 
ciently high strains due to the formation of locally comphant 
regions. In essence, the system exchanges a large elastic strain 
energy per unit length throughout its length to produce a local- 
ized internally disrupted region. The energy cost of creating 
that localized disruption is then compensated by the collapse 
of strain energy into this more compliant part of the filament. 

In this £irticle, we examine this minimal extension of the 
WLC model - the helix/coil wormlike chain (HCWLC) [30]. 
This extension requires us to introduce four material parame- 
ters. These are the bending moduli of the filament in its two 
states, the (free) energy cost per unit length to transform the 
filament from less compliant to more compliant, and a "do- 
main wall" energy associated with extra (free) energy cost of 
a boundary between the two internal configuration states of 
the filament. These latter two parameters are well known in 
the extensive literature on the helix coil model reviewed in 
Ref. [35]. We then apply this model to examine quantitatively 
the recent results of Qu et al. [26] in which they constructed a 
bent piece of dsDNA much shorter than its persistence length 
by attaching the free ends of this oligomer to single-stranded 
DNA (ssDNA). In essence, this construct is akin to a bow as 
used in archery, where the tension in the bowstring (ssDNA) 
is balanced by the bending of the bow (dsDNA); see fig. la. 
By measuring precisely the elastic energy stored in the bent 
dsDNA, Qu et al. obtained a nonlinear torque vs. angle curve 
for the dsDNA oligomer. We use the HCWLC model to in- 
terpret these results. It is important to note that the dsDNA 
in the Qu construct has a nick, a point where the covalently 
bonded DNA backbone is broken. This defect serves to lo- 
calize the generation of the compliant region to the center of 
the construct, a feature that we exploit when developing our 
model. 

The remainder of this article is organized as follows. We 
first develop the general theory in §n, and then examine the 
results of applying a constant force to such a composite fila- 
ment (§111). We then apply the theory to the experiments of 
Qu et al. [26]. We show that one can account for the entire 
measured elastic energy of their construct using well-known 
values for the thermal persistence length of single and double 
stranded DNA by adjusting the two remaining helix coil pa- 
rameters in the model. We find the required fitting parameters 
are consistent with values obtained from other experiments 
suggesting that the HCWLC model provides a consistent min- 
imal description of the nonlinear mechanics of DNA, and pre- 
sumably other filaments having internal structural transitions 
that can be treated as a two-state internal variable. Finally we 
summarize our results and suggest future directions (§IV). 



II. MODEL 



In order to capture the essential characteristics of filaments 
with internal degrees of freedom we use a multi-domain Eu- 
ler elastica with a prescribed energetic cost for melting (see 
definition of terms in fig. 2). Our assumption is that, under 
the influence of external forces or torques, a "bubble" (or soft 
region of some variety) may nucleate at the center of the fila- 
ment, causing a reduction in bending strain at the cost of intro- 
ducing the melted region. We treat the two-domain filament as 
lying in the plane; this will be the lowest energy configuration 
for the Qu construct, which cannot support twist energy due 
to the attachment of the ssDNA. Although there should be out 
of plane fluctuations that contribute to the entropy of the poly- 
mer, these will not efifect the bending of the stiff region, and 
thus is irrelevant to the melting transition under study. In 2D 
we can parameterize the deformation using two tangent an- 
gles: (p{s) and 9(s) for the soft and stifle regions, respectively. 
We let the soft (stiff) region have a bending rigidity denoted 
by /c<(/f>), and thus the bending energy of the system is given 
by 




where 2L is the total length of the filament and 2£ is the length 
of the melted region; we assume here and throughout this 
work that the system is symmetrical about its midpoint. We 
augment this energy functional with the internal free energetic 
cost for melting: 



Ei„, = m + J, (2) 

where the parameter H is the excess free energy per unit length 
to be in the thermodynamically disfavored state. For dsDNA 
(at zero imposed stress), this is the cost of denaturing the fila- 
ment. For a CNT or other tubelike structure is the cost of in- 
troducing ripples or wrinkles. The parameter J is the domain 
wall energy cost of introducing a boundary between the two 
intemal states of the filament. Because we allow for a con- 
tinuous melted fraction but only one domain, the energetic 
cost of J simply shifts calculated criteria for melting. 

Without an external stimulus the rod is trivially straight 
with no melted region, but under the action of an external 
source non-trivial behavior will appear By imposing an exter- 
nal compression force F the total energy functional is modified 
such that 



Em = Ebend + Eint ~ ^ COS (f>ds + j COS 0dsj , (3) 
leading to equations of equilibrium given by 
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FIG. 1 : (Color online) Internal structure in nano filaments can lead to non-trivial mechanical behavior, a) A hybrid loop of dsDNA and ssDNA 
can be schematically represented as a filament composed of a stiff (red, solid) region and a softer (blue, dashed) melted region, constrained 
at it's ends with a spring, b) Under external bending, CNTs develop rippling or kinking modes of deformation that can also be described by 
postulating phase coexistence between two regions (reproduced from [9] with permission). The rippled region has an effectively lower bending 
rigidity, but there is an energy cost associated with inducing the ripples. 
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FIG. 2: (Color online) Schematic of the melting elastica. 
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(4) 
(5) 



These are the Euler-Lagrange equations for the classical prob- 
lem of an elastic filament subjected to a constant longitudinal 
force F. 



dd> do 

ds ds 



(6) 



1 dd)'- 1 de^ 

-f (0 - F COS m-h^ -K^ :r - ^ cos e{{) (7) 

1 ds 1 ds 



m = 

de. 



ds 



(L) = 



(8) 
(9) 



Physically the first two conditions correspond to the continu- 
ity of torque through the interface between soft and stiff re- 
gions, as well as continuity of energy density. From a more 
formal standpoint, these conditions enforce continuity of the 
canonical momentum and Hamiltonian, and are sometimes re- 
ferred to as the Weierstrass-Erdmann comer conditions [36]. 
The third boundary condition comes from the assumed mir- 
ror symmetry of the system and that there can be only one 
region of melting, located on this axis of symmetry. The 
fourth boundary condition is from our assumption that the 
force compresses the filament without imposing a torque on 
the end. General formulations of complicated filaments, in- 
cluding shape-memory alloys [37], adhesive fibers [38], and 
botanical branches [39, 40] often include situations where the 
continuity of the the conjugate momenta and the Hamiltonian 
are augmented by a jump condition of some variety; in addi- 
tion to this the tangent angle need not be continuous through 
the interface in the presence of shear stresses over the cross- 
section. For our present purposes it is sufficient to consider 
that such stresses are contained phenomenologically through 
the domain wall cost J, and we impose (p{l) - 0(() to com- 
plete our formulation. 



The appropriate boundary conditions are found by demand- 
ing that the boundary terms associated with the variational 
derivative vanish: 



By nondimensionalizing lengths by L, and rearranging so 
that the equations of equilibrium are expressed solely in terms 
of non-dimensional parameters we get the following 
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^+AVsin,^ = (10) 
cP-0 

-^+/sin0 = O, 

where A = ^/k^/k^ is a measure of how drastically the two 
regions differ in rigidity, and / - F/F^, with F,. = k^/L^ pro- 
portional to the critical force required for classical Euler buck- 
ling. Similarly, the boundary conditions become, in terms of 
these variables, 



as as 


(11) 


2 ds 2 ds 


(12) 


m = 


(13) 


d0 

:r(i) = o, 

ds 


(14) 



where e - IjL and we have defined h - HL^ Ik> as a di- 
mensionless measure of the energy penalty for introducing the 
melted region. The reduced energy cost for introducing a do- 
main wall is j - JjFcL. Now, given a set of {/i, j. A, /) we 
obtain the size of the melted region that minimizes the energy, 
denoted by e* . 
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FIG. 3: (Color online) Melted fraction for the minimum energy con- 
figuration as a function of applied force. For small loads the filament 
does not melt, but at a critical point (dependent on the melting cost 
h) the optimal value for the melted fraction increases continuously 
from zero. At yet higher values of the force there is a discontinuous 
jump in the stable equilibrium point; these solutions, while interest- 
ing, are of limited value for our present study and their presence has 
been included merely for completeness. Parameter values used are 
A = 8and ; = 0.15. 



Qu et al. 



III. RESULTS 

The material parameter space for this particular model is 
formally fairly expansive, but with an eye towards the DNA 
problem, and armed with the intuition gained from crumpling 
and rippling transitions, we can narrow our search to certain 
regimes. We concentrate on the parameter region A > 1, since 
only if there is a significant bending energy gain upon change 
from melted to unmelted will the nonlinear elasticity regime 
be accessed. For DNA we estimate the persistence lengths to 
be ~ 50nm and ~ O.lnm for dsDNA or ssDNA, respectively; 
as a result A ~ 8. We consider a broad range of h and j, 
although the energetic cost for denaturing dsDNA has been 
reported as * 1 - l.SksT/bp, while the domain wall cost is 
J ~ I - IQksT; these two values correspond to dimensionless 
parameters of h ~ 3 and ; ~ 0.1 [20, 30, 35]. For dsDNA 
we attribute h mainly to the breaking of hydrogen bonds in 
localized base-pair melting while in CNTs this energy scale 
corresponds to the work required to produce the short wave- 
length rippling as shown in fig. lb. 

We present the solutions to two problems. First we consider 
a constant force problem in which a fixed compressive load is 
attached to the ends. In the second we turn to the experimental 
system of Qu et al. To do so we include the mechanics of a 
flexible ssDNA chain connected to the two ends of the much 
stiffer dsDNA filament. We then compute the elastic energy 
of this construct as a function of the degree of polymerization 
of the ssDNA and compare our results to the measurements of 



A. Constant force ensemble 

The simplest result, and perhaps the most intuitively acces- 
sible, stems from examining a constant force ensemble. Un- 
der a constant load, as is generally the case in experiments on 
CNTs, we determine the equilibrium conformations resulting 
from the melting elastica model. One particularly attractive 
reason for examining this formulation is the ease of computa- 
tion, as the equations of equilibrium can be solved directly by 
using Jacobi elliptic functions (see Appendix). 

Increasing the applied load (see fig. 3), we observe a dis- 
continuous melting transition under stress in which the frac- 
tion of the chain in the melted state jumps from zero to some 
finite value at a critical value of the force fc - g(h, j, A). The 
dominant balance between the energy required for melting, 
he and the work done by the applied load, ~ fcosO yields 
the scaling of the critical force to be fc ~ h. At even higher 
values of the applied force new solutions appear in which the 
filament makes one or more loops (see fig. 3), although these 
solutions are unlikely to be experimentally observable. 

We now turn to an analysis of the experimental construct of 
Qu et al. in which the force applied to the dsDNA is provided 
by attaching a flexible ssDNA chain to the ends. We compute 
the experimentally observable elastic energy of the system as 
a function of the arc length of the ssDNA, i.e. the degree of 
polymerization. 



5 



B. Gaussian chain coupled to the HCWLC 

We treat the ssDNA, which is many times longer than its 
thermal persistence length t^, as a Gaussian chain. Includ- 
ing the lowest order nonlinear corrections we write the elastic 
energy contribution as [2] : 



9kBT 



3C 



where 



- I cos (f>ds + I 



cos dds 



(15) 



(16) 



is the projected length of the filament in the horizontal direc- 
tion, Ns is the number of base pairs in the ssDNA. In the limit 
of small stretching we can modify the energy function given 
above by the addition of a spring potential: 



E,„t — E, 



bend 



+ E„ 



.+ -k 
^ 2 '^."-'^11 ' 



(17) 



9kBT 



is the spring constant associated with the har- 
monic part of the stretching energy of the single DNA strand. 
For larger values of Ly, approaching non-negligible fractions 
of the total ssDNA length, including the anharmonic contribu- 
tions of the Marko-Siggia stretching potential would be nec- 
essary. However, we verify a posteriori that we can reproduce 
the nonlinear elastic effects of Qu et al. using the harmonic 
contribution alone. We discuss further the role of elastic non- 
linearities of the ssDNA chain when considering end-to-end 
length fluctuations of the construct in section §111 C. 

Now, instead of a constant force ensemble we must deal 
with the integral constraint on the ends of the filament. Be- 
cause the stretching energy of the ssDNA enters the energy 
functional as an integral equation, it is simpler to minimize the 
energy by constructing the energy landscape numerically and 
minimizing over the parameter space of choice. We nondi- 
mensionalize the energy by /L, such that E/FcL ~ 1 cor- 
responds roughly to lOksT for the relevant values of /<•> and 
L. Ideally we wish to compare to the experiments of Qu et 
al., and thus will examine the elastic energy as a function of 
Ns (i.e. the shorter the strand, the stiffer the spring). In this 
scenario, we expect that the energy cost Hf will correspond 
roughly to the cost for melting a "bubble" of DNA, and thus 
should be on the order of several kgT; similarly, the dimen- 
sionless domain wall cost j = J/fL ~ (9(10"'). The persis- 
tence length of dsDNA is approximately 50 nm, while that of 
ssDNA is ~ Inm, leading to a value of A = ^/k^/k^ ~ 8. 

Before discussing the experimental case specifically, we ex- 
amine more generally the stress induced melting phase dia- 
gram for the Qu constructs using a fixed A = 8, a value appro- 
priate for single- and double-stranded DNA. We also choose 
a domain wall energy j ~ 0.15, well within the range con- 
sidered appropriate for DNA [20, 35]. The phase diagram, 
shown in fig. 4, is spanned by the degree of polymerization 
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FIG. 4: (Color online) Phase diagram for melted elastica constrained 
by a spring. For low values of H the entire filament melts under the 
action of the spring (black region on left of diagram). However, for 
all other values there is a sharp transition between the unmelted and 
partially melted phase. Low values of A', indicate a very stiff spring, 
and thus the filament can remain partially melted for larger H. On the 
other hand, for larger A', the spring acts only weakly, and we see a 
very swift transition from completely melted to completely unmelted 
over a relatively small region of h. The dashed (red) line designates 
the phase boundary derived from scaling arguments. The parameters 
used are A = 8 and j = 0.15. 



of the ssDNA, Ns, and the nondimensionalized free energy 
cost per unit length for being in the more compliant melted 
phase. The contours reflect equilibrium fraction of the melted 
dsDNA with darker colors denoting a larger melted fraction. 
As one might expect, we observe three distinct phases: (i) a 
completely melted filament obtained for small values of h and 
Ns, (ii) a partially melted or phase coexistence region for in- 
termediate values of the two control parameters, and (iii) for 
sufficiently high A^, (such that the ssDNA exerts a low com- 
pressive load on the dsDNA) or high h (such that the free en- 
ergy cost for melting is prohibitively high), the chain does not 
melt at all. We observe a jump in the length of the molten 
region e ~ / at the phase boundary, and a coexistence bound- 
ary that scales as Ns ~ hP; here z and /3 are exponents that 
can be estimated from energy considerations. For the partially 
melted phase the dominant balance in the energy is between 
the melting energy and the spring potential, i.e.. 



he + j 



Ns 



-Ll 



(18) 



We immediately find the exponents z - I and /? = -!, lead- 
ing to a phase boundary scaling of A^j ~ Ijh, with a position 
that shifts depending on the value of j. The dashed (red) line 
in figure 4 shows the scaling behavior is consistent with our 
calculated phase boundary. 

Our primary result focuses on this sharp transition between 
the partially melted and unmelted filament; this generic be- 
havior leads to the ability to recreate the nonlinear elastic re- 
sponse measured in tightly bent dsDNA by the experiments of 
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Qu et al. As seen in figure 5, for short springs the energy be- 
haves hnearly, while the response for longer springs follows 
a power law. Both of these regimes are consistent with the 
data measured by Qu et al. (red dots in 5), and we can use the 
parameters of the model (the melting cost H and domain cost 
j) to fit to these data. We find for dsDNA length 6nm that 
h « 2.1, and the value for j » 0.15. For the case of A^^; = 24, 
the length of the dsDNA is A" « 8nm and we find values of 
h » 0.78 and j « 0.12. Values for h and j reported in the liter- 
ature range from /i ~ 1.5 - 3 and 7 ~ 0.05 - 0.5 [20, 30, 35]. 
Our model captures the possible sequence dependence of the 
melting cost H, as well as displaying a relatively low value for 
the longer dsDNA piece, most likely a result of the destabiliz- 
ing nick at the center of the construct. 



variables to the original axes, Uij^j = ^„ we write 



<m-Lf^f >=kBT 



21 

ri 



u. 



2 ^ 

22 

r2 



(21) 



Taking experimental parameters {Ns = 15) and our fitted value 
oih = 2.1 {Nd = 18) we predict that the end-to-end fluctua- 
tions of the Qu construct wiU be quite large. In fact we find 
that 



^<(L||-L[l''>)2> 



r(0) 



0.6 



(22) 



C. Effect of fluctuations on melted region 

We expect the phase coexistence region to be characterized 
by an anomalously soft elasticity. If the end-to-end force on 
the chain is increased, the chain can accommodate by a com- 
bination of bending and increasing the size of the melted re- 
gion. This soft elasticity associated with the coexistence re- 
gion should then be reflected in the equilibrium statistics of 
end-to-end distance fluctuations. In this section, we use our 
HCWLC model with a melted region to probe these fluctua- 
tions and make predictions yet to be tested. 

To be specific, we compute the end-to-end fluctuations of 
a Qu construct by examining the Gaussian fluctuations of the 
end-to-end distance of the dsDNA chain for a fixed polymer- 
ization index A^, of the ssDNA, which we continue to treat 
as a linear spring. Expanding the energy of the system about 
the point of mechanical equilibrium (€^^\ ^y'^')' write the 
energy of the system as 



1 /)2<^y(/) 

2 d^id^j 



(19) 



where ^1 - Ly - and ^2 - e - e^"' measure deviations 
of the size of the melted region and end-to-end distance, re- 
spectively. The entire energy landscape of the Qu construct is 
computed numerically and shown in fig. 6. The quadratic ap- 
proximation to the energy is valid near the minimum (shown 
as the yellow filled circle in fig. 6) and can be simplified by 
rotating to the principal axis frame; due to the complex elastic 
nonUnearities of the system, it is more convenient to compute 
both the principal axes £,\ , ^2 and corresponding curvatures nu- 
merically. In the principal axis frame the Hamiltonian takes 
the form 



[1 .2 

2 ^1 + 2 ^2. 



(20) 



where ri,r2 are the principal curvatures associated with 
the energy landscape. A straightforward application of the 
equipartition theorem yields the mean squared fluctuations of 
^1,^2- Using the rotation matrix U relating the principal axis 



For an unmelted dsDNA chain in the Qu construct, the prin- 
cipal curvatures are significantly higher, restricting the end-to- 
end fluctuations to a much smaller amount: the classical result 

for the mean square end-to-end distance < >= I? fr,{LI(p), 

wifli /z)(x) = {\-x + e '')/x^ flie Debye function [41, 42]. For 
the dsDNA in the Qu construct L<^ {p and thus the root mean 

square end-to-end fluctuations -y/<(L^-i^i|)^~>/i' ~ g^/^p ~ 
0.02. 

One can understand the dramatic softening of the system 
in the partially melted state as follows. In the coexistence 
region, changes in the total bending energy can be partially 
compensated by adjusting the boundary between the compli- 
ant melted region and the stiffer dsDNA. The analogy to the 
compressibility of a mixture of gas and liquid in the two-phase 
coexistence region is apt. In both cases the incremental work 
associated with an increase in the mechanical load remains 
fixed (instead of increasing), and that work converts material 
in one phase to the other. 

The definitive test of the prediction of the enhancement of 
the end-to-end distance fluctuations in the Qu construct due 
to local dsDNA melting is to compare the width of the end- 
to-end distance distribution for Qu constructs with various ss- 
DNA lengths, A'^j. Based on the above analysis, we expect one 
to observe a jump in the width of that distribution as A^, is 
decreased past A'*, where A'* is the critical length of ssDNA 
below which local melting first occurs. This length can be in- 
dependently extracted from the data by finding the break in 
the slope of the £tot vs. A^^ curves, as shown in Fig. 5. 

Since the separation of the dsDNA ends in the Qu con- 
struct is less than lOnm, FRET pairs attached to the ends of 
the dsDNA should be accurate reporters of the needed end-to- 
end distance [43, 44], and there are well understood methods 
available for attaching FRET donor acceptor pairs to DNA. 
By using single pair fluorescence resonance energy transfer 
(spFRET) [45], one should be able to extract data on the distri- 
bution of this end-to-end distance in freely diffusing Qu con- 
structs in solution. In fact, this technique has been success- 
fully used to observe the unfolding of a DNA hairpin [45], 
and thus may be extended to observe the end-to-end distance 
distribution of the dsDNA in the more complex Qu construct. 

It is important to note that our calculation makes two sim- 
plifying assumptions. The first is that the ssDNA can still be 
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FIG. 5: (Color online) Total elastic energy of the ssDNA spring and dsDNA system, as a function of the number of ssDNA base pairs A^j with 
the number of base pairs A'^ = 18 shown in part (a) and A'^ = 24 shown in part (b). The model prediction is shown as the solid (blue) line. The 
filled (red) points are data from [26], fit using the parameters A = 8, j = 0.l5,h - 2.1, for A^^ = 18 and A = 8, /j = 0.78, ; = 0.12 for A'^ = 24. 
In the low A', region the spring is relatively stiff, leading to a higher melted fraction in the dsDNA and an efi'ective strain softening of the entire 
construct. 



treated as a linear spring. It appears that this cannot be valid 
of the full predicted range of extensions, as the fluctuations in 
L|l can become large. The experimental system may in fact 
be somewhat stifFer for this reason, but the significant fluc- 
tuation enhancement of the two-phase regime over that of an 
unmelted Qu construct remains valid. Nonlinear corrections 
to the ssDNA energy are not necessary to understand the equi- 
librium state of the experiments and we do not consider their 
effect on the energy landscape here, although such corrections 
are straightforward to pursue numerically. The second as- 
sumption used here is that the internal state variable remains 
in thermal equilibrium throughout the fluctuations. We expect 
that this is so, but it may be that the melting/unmelting of the 
dsDNA will be sequence-dependent and that some sequences 
may have slower kinetics. 




0.85 0.7 0.55 0.3 0.15 0.0-0.15 
L 
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IV. DISCUSSION AND SUMMARY 

The Qu construct presents an important experimental probe 
of high energy (small radius of curvature) DNA bending me- 
chanics. In such a regime it is now clear that DNA bend- 
ing becomes nonlinear and thus cannot be fully understood in 
terms of the WLC. This nonlinear softening of dsDNA due to 
accessing extra local degrees of freedom in high elastic en- 
ergy deformations, which we associate with chain melting, is 
not in itself surprising. Nevertheless, understanding how to 
account for such effects in an analytically tractable model is 
essential for the further study of the biophysics of highly de- 
formed dsDNA or DNA under large external loads. To this 
end, the HCWLC model and it variants have been proposed 
as natural extensions of the WLC model applicable to such 
higher energy mechanics of DNA and other filaments with in- 
ternal degrees of freedom, including alpha-helical polypep- 
tides [20, 30-32, 34]. 

We have shown here that a restricted elastica calculation 
based on the HCWLC is able to account in quantitative detail 
for the observed nonlinear elasticity of dsDNA as measured 
by Qu and collaborators. To fit the data from the Qu construct 



FIG. 6: (Color online) Energy contours in the space of melted region 
and projected length. The (yellow) filled point represents the zero 
temperature energetic minimum, the (white) dashed lines show the 
directions of principal curvature, and the (red) elliptical contour cor- 
responds to kgT energy fluctuations away from the minimum, calcu- 
lated by expanding the energy landscape to quadratic order about the 
energy minimum. While the melted region remains relatively con- 
stant, the applied force undergoes large fluctuations in the presence 
of Gaussian thermal excitation. 



using dsDNA with Nd - or A^^ - 24 base pairs, we fix the 
thermal persistence lengths of the dsDNA and molten DNA to 
be 50nm and ~ Q.lnm consistent with measurements of dou- 
ble stranded and single stranded DNA respectively [46, 47]. 
Using those fixed elastic constants, we find the data of Qu 
and collaborators can be fit by choosing H » .96ksT/bp (for 
Nd - 18) or H X Q.lkgT /bp, on the same order with that 
predicted and measured elsewhere [20, 35]. The effect of the 
domain wall cost J, found to be 7 3kBT, is relatively in- 
dependent of our determination of H and compares well to 
previous helix/coil values of between 1 - lOk^T determined 
from temperature induced dsDNA melting studies discussed 
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in [20, 30, 35]. For the parameters used, the energy cost as- 
sociated the bubble dominated by H term , i.e., H{ > J. For 
systems in which the compliant domain has an even softer 
bending modulus, we expect the arc length of the molten re- 
gion to be smaller. In this case, the domain wall energy may 
dominate the length-dependent contribution to the free energy 
cost of the bubble. 

In addition to showing that the present elastica calculations 
can account for the dsDNA bending elasticity as measured 
by Qu and collaborators, we suggest that the presence of a 
molten bubble leads to anomalously large end-to-end distance 
fluctuations of the Qu construct in the "two-phase" regime as- 
sociated with the experimental range Ng < 25 - see Fig. 5. 
These large fluctuations result from the correlation between 
the bending of the dsDNA and change in the arc length of the 
molten bubble. We propose the observation of a qualitative 
change in the end-to-end fluctuation spectrum of the Qu con- 
struct at the kink in the energy versus A^, curve represents a 
key additional test of the theory. The calculations of section 
IIIC provide a quantitative prediction of that effect. In those 
equilibrium calculations, we have assumed that the internal 
degrees of freedom of the dsDNA remain in thermal equilib- 
rium on the time scale of the observed end-to-end fluctuations. 
While we expect this is reasonable at present, one can address 
the potential for nonequilibrium kinetics in the melting and 
"refreezing" of the dsDNA in a manner analogous to that use 
in the context of elastic deformations of shape-memory al- 
loys [37]. 

The principal restrictions made to the full HCWLC model 
in this calculation were that the apphed load produces at most 
one molten region of the dsDNA and that this unstructured 
part of the molecule is nucleated at the midpoint between 
the ends of the dsDNA; it also grows/shrinks symmetrically 
about that point in response to changes in the applied load. 
In the case of the experiments in question, these assumptions 
are reasonable. In the case of the former, the domain wall 
energy associated with producing extra molten regions 2J is 
sufiiciently high so as to strongly suppress such fluctuations. 
To justify the latter, one must recall that the Qu construct ds- 
DNA is nicked at the midpoint, providing a nucleation center 
for melting. Moreover, the dsDNA is sufliciently short that 
the entropic contribution to the free energy associated with 
the translational degree of freedom of the molten bubble is a 
subdominant correction, as discussed in Ref. [30]. 

We also do not consider the effect of quenched disorder in 
this set of calculations. In DNA one expects the local melting 
temperature to be suppressed in regions where there is a high 
density of AT base pairs. Moreover, there is some evidence 
that the bending modulus of the molecule also depends on se- 
quence. We have not explored this issue, but one may imagine 
that the position of the kink in Fig. 5, corresponding to the nu- 
cleation of a molten region and determining our fit for H, may 
depend in detail on the sequence. In addition, the application 
of this theory to microtubules will have to eventually confront 
the role of quenched disorder since these structures are known 
to possess a number of defects in local protofilament structure 
and these are believed to influence both the zero stress config- 
uration of the filaments and their local elasticity [48-52]. The 



same can be said for carbon nanotubes [4, 9, 53-55]. 

This model, although focused on the bending of biopoly- 
mers with complex internal structure, can be applied to other 
types of soft matter elasticity. Our variational formulation 
is similar in many respects to that employed to describe 
the peeling apart or adhesion of two or more elastic fil- 
aments [38, 56, 57] and to discuss elastocapillary coales- 
cence [58]. As such, it is possible that further application of a 
multiple domain filament with a measurable cost for peeling 
could be used in such systems to simplify the analysis. 
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Appendix A: Elastica solutions in terms of elliptic fmictions 

The elastic energy functional for an elastic rod under uni- 
form compression is 



Jo 



K idey 



ds. 



(Al) 



where L is the length of the rod, k is the flexural rigidity, 6 
is an angle that parametrizes the rods tangent vector, 5 is an 
arc length coordinate, and F is the compression force. A first 
variation of this integral yields the Euler-Lagrange equations 
for the so-called "elastica": 



d 6 

K—r+FsiQe = 
ds^ 



(A2) 



This equation, despite being nonlinear, admits a first integral 
in 6, indicating that energy is conserved: 



1 Ide 

2\ds 

2 



- F cos 6 = -F cos a, 



(A3) 
(A4) 



where we have chosen the constant of integration such that 
the angle at s = L is identically equal to a, for rods with an 
inflection (dO/ds = 0), and the modulus k is chosen so that 
filaments without inflections will have a tidy formulation (see 
below). There are two main classes of solution, referred to as 
inflectional (I) and non-inflectional (NI); the former involve 
solutions that have dO/ds = somewhere on the filament, 
while the latter do not [59]. 

Scaling lengths by the filament length L, the solutions for 
the curvature dOjds everywhere on the filament are given in 
terms of the elliptic functions cn{s) and dn{s) [60]: 
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^ =2A:Vfcn ( Vf 5 + /?•), (7) 



(NI) 



(A5) 
(A6) 



Here we have defined the reduced force / - FL^/k, the com- 
plete elUptic integral of the first kind K, and the Jacobi elliptic 
functions are defined in the following manner: 



,k)= r — 

Jo Vi 



dt 



■ sin^ t 



cn{s) = cos ^ 



(A7) 

(A8) 
(A9) 
(AlO) 

(All) 



The angle and the Cartesian coordinates are now given by 
various integrals, 



Jo 
Jo 



de(s') 



ds' 



ds' 
cos e(s')ds' 



yis) 



f sin 9is')ds', 
Jo 



(A12) 
(A13) 
(A14) 



all of which can be expressed in terms of Jacobi elliptic func- 
tions or elUptic integrals of various kinds. 

Using these general solutions, we need only apply bound- 
ary conditions in order to find particular shapes of an elastica 
under compression. For the case of the melting elastica under 
a constant compression force we can define the inner melted 
region as a non-inflectional solution, and then match to the 
outer filament, which satisfies the inflectional solution. 

The boundary conditions are: 



4>(L) = 

de 

-(0) = 
ds 

do dd> 
ds ds 



(A15) 
(A16) 

(A17) 



These three conditions come from mirror symmetry at s = 
L, torquelessness at s = 0, and torque continuity where the 
regions change from soft to stiff (s = £). For the time being, 
we assume that the softer region has a fixed length {, so that 
our problem reduces to matching two separate solutions for 
the elastica: <p(s) satisfies the conditions for (NI) solutions, 
while 6(s) corresponds to (1) boundary conditions. 

The inflectional solution for dO{s)lds is given by 



J-^='^ki^ cni ^s + Kiki)), 0<s<€ (A18) 

Similarly, the non-inflectional solutions for d<pis)/ds are 
given by 



e< s <l 



(A19) 



where e = ([L is the dimensionless position of the domain 
wall. Given these two solutions we now need only find the 
value of the elUptic moduli ki and k2, which can be found by 
applying the matching condition (in dimensionless form): 



(e)=^(e), 
as ds 



and the condition of force continuity: 



(A20) 



i as 



(e) 



2H. 



(A21) 



Using equations A18 and A19 the matching condition can 
be formulated as a transcendental equation for ki and k2: 



Akik2 cn [fe + K{k\)) = -dn | 



A/(l - e) 



{All) 



The relationship between the two elliptic moduh is given by 
the second boundary condition, and thus given values for /, 
e, A and H we can now find solutions by solving two tran- 
scendental equations. The minimum energy solution can then 
be found by minimizing the energy with respect to the melted 
region e. 
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